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A STABLE AND CONSERVATIVE INTERFACE TREATMENT OF ARBITRARY 

SPATIAL ACCURACY* 

MARK H. CARPENTER*, JAN NORDSTROM*, AND DAVID GOTTLIEBS 


Abstract. Stable and accurate interface conditions arc derived for the linear advection-diffusion equa- 
tion. The conditions are functionally independent of the spatial order of accuracy and rely only on the form of 
the discrete operator. We focus on high-order finite- difference operators that satisfy the summation-by-parts 
(SBP) property. We prove that stability is a natural consequence of the SBP operators used in conjunction 
with the new boundary conditions. In addition, we show that the interface treatments are conservative. 

New finite-difference operators of spatial accuracy up to sixth order arc constructed: these operators 
satisfy the SBP property. Finite-difference operators are shown to admit design accuracy (/'/'‘-order global 
accuracy) when (p - l^-order stencil closures are used near the boundaries if the physical boundary con- 
ditions arc implemented to at least p th - order accuracy. Stability and accuracy arc demonstrated on the 
nonlinear Burgers equation for an twelve-subdomain problem with randomly distributed interfaces. 

Key words, high-order finite-difference, numerical stability, interface conditions, summation-by-parts 

Subject classification. Applied and Numerical Mathematics 

1. Introduction. Higher order and spectral schemes arc ideally suited for resolving problems for which 
high resolution is essential. Computational acroacoustics (CAA) and computational electromagnetics (CEM) 
arc two such fields that require high accuracy to resolve the vastly disparate length and time scales involved. 
High-order (spectral) schemes easily outperform low-order schemes on simple problems in which the physical 
domain is smoothly mapped onto the computational space. The spatial convergence rates of these schemes 
allow satisfactory results on relatively coarse grids. 

At least two fundamental obstacles presently limit the use of high-order schemes. The first one is the 
lack of nonlinear robustness exhibited by high-order formulations. Under resolved features in the solution 
and inappropriate numerical and physical boundary conditions are the primary causes. A second limitation 
is the difficulty in applying high-order formulations to complex geometries. Often, the generation of a grid 
around a complex configuration is the most difficult aspect of the solution procedure. Further constraint 
of the grids so that they are smooth to higher order (necessary to attain design accuracy for high-order 
methods) severely complicates grid generation around complex configurations. 

Many high-order practitioners advocate a fully unstructured approach to grid generation. This approach 
simplifies the grid-generation procedure considerably for complex configurations. Finite-clement techniques 
arc an example of the fully unstructured schemes that are routinely used on complex geometries. An 
alternative to fully unstructured methods is the semistructured approach, in which the solution domain 
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is broken into the union of piecewise smooth subdomains. Each subdomain is discretized with a stable 
formulation, and the resulting multiple domains are patched together globally. This technique has been 
successfully used by Kopriva [1], and more recently by Hesthaven and Gottlieb [2]. 

The approach for designing interface conditions developed in this work is equally valid for the unstruc- 
tured and semi-structured approaches in multiple spatial dimensions. The interface conditions arc determined 
entirely by accurate left and right state data along the interface, and do not depend on the source of the 
data. For simplicity, however, we focus on the interface matching conditions necessary to maintain stability 
and accuracy in one spatial dimension. We demonstrate the technique for both spectral and high-order 
formulations. 

In section 2, we define and describe semidiscrete operators that satisfy the SBP convention. In section 
3, we introduce new interface boundary conditions for multiple domains. In section 4, we show that the new 
conditions are conservative across interfaces. In section 5, we consider specific examples of the stability and 
accuracy of finite-difference schemes. In section 6, we present the conclusions. Finally, in the appendix we 
present the stencils used for fourth- and sixth-order finite-difference schemes. 

2. Spatial Discretizations. The stable interface conditions presented in this work are valid for spa- 
tial discretizations of arbitrary accuracy. To achieve this generality, the spatial discretizations must be of 
a specific form. Fortunately, most numerical schemes can be put into the required form with only mi- 
nor modifications. To be more precise we consider discrete spatial derivative operators with the following 

properties: 

2.1. First-derivative properties. 

1. The first derivative operator defining the numerical derivative u x = [(^)oi”^(^)n] is 

Pu x — Qu = 0 

(2.1) Fv x - Qv = PT e , 

where u = [u 0 {t),ui{t), ...,u N {t)} T , v = [t>(x 0 , t), v(:en, t)] T and v x = [(gf )o, (||)n] T (The 

vector v is the exact solution.) The truncation error T e satisfies |T e | = 0(Ax) m where the quantity 
Ax is defined as the maximum distance between any two neighboring grid points. 

2. The matrix P is symmetric and positive definite (A x)p I < P < (Ax)q I where p and q arc 
independent of N with p > 0 and q > 0. 

3. The matrix Q is nearly skew symmetric and satisfies the property Q + Q T = D, where the diagonal 
matrix D has the form d M = [-1, 0, ..., 0, 1] for i = 0,1,. ..TV. Furthermore, Qo,0 = “2 and 

Qn<n — 

A spatial operator in the form of equation (2.1), which satisfies properties 1 through 3, is referred to as 
an SBP operator [3]. All SBP operators automatically lead to an energy estimate for periodic solutions to 
the linear advection-diffusion equation. In the finite-domain case, an energy estimate exists when an SBP 
operator is combined with specific boundary treatments. 

Discretization operators that satisfy the SBP framework are remarkably general. Kreiss and Scherer [3] 
first suggested the use of SBP spatial operators in the context of second-order central-difference schemes. 
In Olsson [4] [5) [6] and Strand [7], high order finite difference operators arc constructed based on spatial 
operators of SBP type. These resulting schemes are strictly stable which means that the growth rate of the 
analytic and semi- discrete solution is identical. 

The precise properties of the matrices P and Q provide a constructive means of formulating boundary 
closures. A discretization begins with a parameterization of several points near the boundary of the required 
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accuracy. The parameters arc then adjusted so that they match the precise requirements of the P and Q 
matrices. Strand [8] used the SBP approach to construct stable fourth- and sixth-order central-differencing 
schemes with boundary closures of the appropriate order. Carpenter, Gottlieb, and Abarbancl [9] extended 
the SBP formalism to compact implicit operators (fourth-order Pade’ operators); Carpenter and Gottlieb [10] 
showed that spectral formulations (Galerkin and collocation) can be cased in the SBP framework. Finally, 
Carpenter and Otto [11] showed that the SBP schemes have a natural interface property, and they used this 
property to derive a class of multiple-domain schemes referred to as “cyclo-difference” schemes. (The earlier 

work [11] required strong imposition of interface data, whereas the present formulation requires only weak 
imposition). 

The SBP schemes naturally arise with centered approximations, for which the spatial operator is skew 
symmetric. A more general class of schemes could be formulated in the form 

(2 ' 2 ) ~ = P~ 1 (Q + T)u 

where the matrix T is symmetric negative definite. The general formulation includes the entire class of 
central and upwind schemes. The upwind schemes arc automatically stable and accurate because they arc 
obtained by adding a symmetric high-order diffusion operator to a stable and accurate SBP formulation. We 
focus, therefore, on the original SBP definition which includes central, compact, and spectral formulations. 

An approach similar to that used on the first-derivative operator can be used for the second- derivative 
operator. For example, one can seek two positive-definite matrices L and R such that 

vxx - L-'Rv = 0(Ax) m 

An obvious choice is to take L = P and R = QP~ l Q so that the second-derivative operator is obtained by 
repeated differentiation with the first- derivative operator. For spectral discretizations, this process of differ- 
entiation is a natural consequence of the polynomial- based discretization technique. This same assumption 
for finite-difference techniques is acceptable but less desirable than other, more compact formulations. A 
second derivative formed from two first-derivative operators is unnecessarily wide and inaccurate and can 

lead to odd-even mode decoupling. For this reason, we seek a second-derivative operator with the following 
properties: 

2.2. Second-derivative properties. 

1. The second-derivative operator that defines u IX is 

( S T M -f- D)Su = 0 

( 2 - 3 ) Pv xx - (~S T M + D)S\ = PT e2 . 

where the diagonal matrix D has the form d t)l = [-1, 0, ..., 0, 1], i = 0, 1, ...N 

2. The matrix M is positive definite: (Ax) ml <M< (Ax) n I , where m and n arc independent of 
N with m > 0 and n > 0. 

3. The matrix S is of the form 

5 0,2 so , 3 

0 

1 0 


0 1 0 

0 1 0 

S N,N - 3 SN,N - 2 S;v,./V-1 


(2.4) 


5= 


(Ax) 


s 0,0 So,l 

0 1 
0 
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where 


Su\o = + O(Ax) 

(2.5) Su\ N = v x {x N ) + 0{ Ax) r 

The matrix 5 is the identity matrix (scaled by the grid spacing) where a discrete representation of 
the first derivative replaces the first and last rows. 

4. The matrix P is that used in the first-derivative operator. 

Explicit forms of the matrices 5 and M are given in the appendix for a second-order explicit discretiza- 
tion. In addition, the matrix S is presented up to sixth order. 

3. Interface Boundary Conditions for Multiple Domains. Consider the Unear advection-diffusion 


equation 

(3.1) 


dU dU_ 
dt ° dx 


= e 


d 2 U 

dx 2 ' 


<1, t > 0. 


Suppose that the equation is discretized by a multi-domain technique such that the interval is divided 
arbitrarily into two subintervals -1 < x < a* and x t < x < 1. On each subinterval, a discretization is used 
that satisfies the SBP properties 1 through 3. We propose implementing the interface boundary conditions 
by using a penalty treatment of the form 


(3.2) 


P t u t + aQiu = eRiu + o\e\ i (u\ x = Xi — v| x =xj + cr 2 eeu[(D/u)| x=Xi (D r v)\ x - Xi ] 

j P r Vt aQ r \ = ei? r V + 0’3<3ri(v|x=Xi ^ix=Xi) "P |x=Xi (^/u)| x== Xtl 


where u is a vector of length M u = [tio(i), Ui(t), .... u M (t)] T defined in the left domain at the points 
XL = [ Xo = = ^)] T and e lt = [0, ...,0, l] r is of dimension M. In the right domain, v = 

[n 0 (t),uj(t),...,^(t)] T is defined at the points x R = [x 0 = Xi,xi, = 1] T and e ri = [1,0, ...,0] T is of 

dimension N. 

The second-derivative matrices P^Rl and Pp'Rr, as well as the first-derivative matrices P t ~ l Qi and 
p-iQ rj a re defined as in section 2. The matrices D; and D r are any operators that approximate the first 
derivative to 0( Ax)" 1 . The obvious first choice would be to use P l ~ 1 Qi and Pp 1 Q r , but this choice is not 
essential for accuracy or stability. (In equation (3.2) we have ignored the physical boundary conditions at 
x = 1 and x ~ -1 for the sake of simplicity. ) 

THEOREM 3.1. Consider the scheme (3.2) for the advection-diffusion equation (3.1). If the matrices 
p h Qi,P r , Q r . Rl and Rr satisfy the first and second derivative properties of section 2 and 


(3.3) 


<^3 — <J\ 


a, a"4 = + 1, 


a 

ai - 2 


€[ 4^ + ^Z ] 


then (3.2) is stable. 

In the proof which follows, we have without loss of generality considered only the interface terms, and 
ignored the terms that arise at the physical boundaries. We assume that the physical boundary conditions 
are implemented by stable and accurate numerical procedures. (Sec Hesthaven and Gottlieb [2] for a possible 
implementation) . 

PROOF: The proof is based on a simple energy estimate. By premultiplying the equations by the vectors 
u T and v T , respectively, and adding we obtain 

+ IMIpJ = 2u T {eRi-aQi)u + 2v T (eRr-aQ r )v 
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+ 2 OiUi(Ui - Vi) + 2 €<T 2 Ui[{Diu)i - (D r v)i\ 
-f 2 a 3 Vi(vi - Ui) + 2ea 4 Vi[(D r v)i - (Diu)i\ 


where ||ix||^ = u r P<u, and we have defined u*, v*, (Diu) v and (D i v) i as u\, 
(D r v)\ x=Xi , respectively. The second-derivative properties of section 2 lead to 


(Au)| I=Ii , and 


( 3 - 4 ) u T R[U < -ai(Diu)i + Ui{Diu)i 

( 3 * 3 ) v T RrV < -a r (D r v)^ - Vi(D r v)i 


where the constants a/ and a r are positive. 

By using the first- derivative properties of section 2 and equations (3.4) and (3.5) and neglecting the 
physical boundary terms leads to 


( 3,e ) p t + IMIpJ < 

where w 2 = [m, vu (A«)j, (A-v)»]> and the boundary matrix B defined by 

(-a + 2cri) — (o*i+cr 3 ) e(l+er 2 ) -ea 2 

— (^iH-<t 3 ) a -1-2(73 — ea 4 e(— 1 + cr 4 ) 

e(l-f-<7 2 ) —ea 4 —2 ea x 0 

-ea 2 e(-l + <r 4 ) 0 — 2ea r 

Straightforward (though tedious) algebra shows that conditions (3.3) yield a non-positive definite matrix Z?, 

thus proving stability. Details are presented in Appendix I. 

In practice, the values of a\ through a 4 are determined as follows. The parameters a r and ct X are 

functions from the numerical method and the chosen grid. The viscous contribution in the constraint 
2 2 

equation o\ < | — e[^ + ^-] is minimized for <j 2 = , yielding the expression < 7 \ < | — e[ 4 ^ a 

The value o\ determines the dissipation at the interface, and also influences the effective CFL of the numerical 
scheme. Values of a\ in the range — 1 < <T\ < 0 provide a compromise between adequate levels of dissipation, 
and acceptable numerical efficiency. 

We have shown that the linking of two domains at an interface with the interface conditions prescribed 
in Theorem 3.1 is stable in a semidiscrete sense for specific values of the penalty parameters <j\ through o 4 . 
The basic methodology can be extended to an arbitrary number of subdomains without complication. The 
only constraint is that the numerical method must satisfy the SBP framework. The methodology does not 
rely on subdomain size and does not require the same SBP operator to be used in each domain. In principle, 
a finite-difference operator of any order, as well as spectral operators on subdomains of arbitrary size, can 
be linked together in a stable manner. Practical details on how to chose o\ through a 4 are included in the 
results section (Section 6). 

In section 2, we presented the general form of second-derivative operators appropriate for this work. We 
then noted two specific derivative operators that satisfy this form. We now show that both choices for the 
matrices R x (and R r ) suggested in section 2 satisfy conditions (3.4) and (3.5) of Theorem 3.1. We start with 
the first option (i.e. R t = Q i P l l Q l ). In this case, the first derivative matrix in (3.2) is D { = P l l Q i . Thus, 
the quantity u T R x u becomes 



^QiP^Qm = u 

= -(pr 1 Qiu) T PiiPr'Qiu) + ui(pr l Qiu)i 
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where we have used the SBP property Q + Q T = D , and have ignored the physical boundary contribution. 
We recall now that Pi > (A x)pi so that 

u t Riu = u T QiP^ l Qiu 

< -(Ax)pi|(J3ju)| 2 + Ui(Diu)i 

Thus, (3.4) is satisfied with = (A x)p L . A similar result holds for Rr with a r = (A x)p r . 

The second choice presented in section 2 for the second- derivative operator P~ 1 Ri is of the form of 
equation (2.3): 

P l Ri = P-\-S t M -f D)S 

For the purpose of proving stability, we relate the two matrices D{ = S. (In actuality, only the first and last 
rows satisfy Di = S. They are, however, the only portions of the matrices that enter the proof.) 

u T Rm = —(Su) T MSu + Ui(Su)i 
< —(Ax) m|Su| 2 + Ui(Su)i. 


Thus, (3.4) is satisfied with cq = (Ax) m. 

4. Conservation at the Interface. The Lax-Wendroff theorem [12] addresses the complexities en- 
countered in solving nonlinear conservation laws. The theorem states that a convergent numerical approxi- 
mation Ui(x , £), computed with a consistent and conservative method, converges to a weak solution of the 
conservation law. Note that discrete conservation is necessary to satisfy the conditions of the theorem. 

A heuristic definition of conservation (commonly encountered by practitioners) describes how the numer- 
ical flux function “telescopes” across a domain to the boundaries. The total quantity of a conserved variable 
in any region changes only as a result of the flux through the boundaries of the region. We, however, rely on a 
broader definition of conservation motivated by the original proof of the Lax-Wendroff theorem. We demand 
that the numerical flux telescope across the domain, and that all moments of the flux against an arbitrary 
test function telescope across the domain. This additional constraint demands an equivalence between the 
weak forms of the continuous and discrete operators. 

We begin by discussing conservation in a single domain. Consider the nonlinear equation Ut + F x = 0 
on -1 < x <1 and t > 0. Note that in the linear case F = all and we obtain (3.1) with 6 = 0. To 
obtain the weak form of this equation we multiply by an arbitrary test function 0(x,i) that vanishes on the 
boundaries. By integrating with respect to space and time we obtain an integral statement of the original 
differential equation: 

f iUixt, - s:l {U<f> t + F<j> x )dxdr = 0 

Now consider the semidiscrete equation given by PU t + QF = 0. Here, we have replaced the spatial 
derivative F x in the continuous case with an SBP derivative operator of order (Ax) r . By multiplying by 
the discrete vector 4>{xj) = 4> T (the discrete analog of integration) and integrating with respect to time, we 
obtain 

<t> T PU\ l 0 - [\u T P<f> t + F T Q<j))dr = 0 

Jo 

Thus, the semi- discrete operator satisfies a weak form similar to that of the continuous operator, and asymp- 
totically approaches the continuous operator in the limit of infinite spatial resolution. The special form of the 
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P and Q matrices present in the SBP operators enables the semidiscrete operator to mimic the conservation 
property of the continuous operator. 

The equivalence between the continuous and semi- discrete operators is more more complicated for mul- 
tiple domains. The conservation property of the SBP operator does not necessarily apply at an interface 
boundary. Under very mild restrictions, however, the SBP interface operators telescope out to the physical 
boundaries, as does the continuous operator. Because conservation is only necessary for the advection terms 
in the advection- diffusion equation, we set e — 0 (see equation (3.1)) and prove conservation for a two- domain 
discretization. We prove conservation for a general nonlinear flux. Note that the penalty parameters for this 
nonlinear case are designated &i and <r 3 . The resulting conservation condition obtained in the nonlinear case 
is slightly different from that obtained in the linear analysis. This difference results from different scalings 
of the penalty parameters. 

THEOREM 4.1. Assume the nonlinear equation = o is valid on the interval -1 < x < 

1 it >0, divided arbitrarily into two subintervals -1 < x < x t and < x < 1. On each subinterval 
a discretization is used that satisfies the SDP framework , and boundary conditions are imposed via penalties 
in the form 


u t + P~ l Q,F(u) = ^ 1 P i - 1 e li [F(u(x i )) - F(v(n))] 

( 41 ) v < + P r _ 1 g r F(v) = - F(u(ii))] 

where u= [wo(t), tti(£), UA/(t)]™ is defined in the left domain at the points xl = [xo = — l,Xi = x,] T . 

and ei, = [0, ...,0,1] T is of dimension M, with similar definitions on the right domain. The discretization 
is conservative provided that the stability condition d 3 = d\ - 1 is satisfied. 

PROOF: For multiple domains, we proceed as shown previously in the single-domain case. Multiplying 
equations (4.1) by the vectors 4> T Pi and 0 T P r , respectively, yields the set of equations 

4> T Pint + <f> 7 QiF(u) = di <t>(xi)(F(u(xi)) - F(i>(x t ))) 

0 T P r v t + <£ T Q r F(v) = &3<t>(xi){F{v( Xi )) - F(u(xi))) 

Using the properties of Qt and Q r we get 


4> T - F T Qi<j>+ 4>{xi)F{u{xi)) — ai(j>(xi)(F{u{xi)) - F(v(xi))) 
<p T PrVt - F T Qr4> - <t>(xi)F(v{Xi)) = a^<f>( Xi )(F{v(xi)) - F(u(xi))) 


By integrating with respect to time and making use of the fact that <f> is continuous at the interface, we get 

<£ r P ; i4 + <A T P r v|‘ = [ {u T P l <f> t + F T Q l <f>)dr 
Jo 

+ f {v T P r <p t -f F T Q r <f>)dr 
J o 

+ [ 4> x F(u{xi)){a i - <r 3 - 1 )dr 
Jo 

A- I <p t F{v(xi)){crz - + 1 )dr 

Jo 


Obviously, the condition <j 3 = d 1 - 1 eliminates the interface terms from the expression and leaves the 
desired weak form of the semidiscrete equation. Thus, the theorem is proved. 
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5. Accuracy of Boundary Conditions. A significant obstacle in dealing with high-order finite- 
difference schemes is the formulation of stable stencils near the boundaries. A uniformly high-order ap- 
proximation should be maintained if possible up to the boundary. In most high-order formulations, ensuring 
uniform accuracy up to the boundaries is difficult when numerical stability must be maintained. Fortunately, 
Gustafsson [13] showed that difference approximations to mixed hyperbolic parabolic equations admit global 
design accuracy when a finite number of points (independent of N) are closed with stencils that are less accu- 
rate by 1 order. For example, a fourth-order interior discretization will asymptotically recover fourth-order 
Z /2 accuracy with third-order closures near the boundaries. 

In this section, we confirm that the physical boundary conditions must be imposed with at least the 
design accuracy in the context of interface boundary conditions. We begin by inspecting equation (3.2) and 
by defining the truncation error as that error committed by substituting the exact solution into the scheme. 
Denote by V\ and V r the projection of the exact solution in the two domains. Substituting the exact solution 
into the first equation in (3.2) yields 

P/Ti c = Pi ^ 4- aQ\V\ — eR[\ \ + — V ri = x .) + cr2&\i{{DiV\)\ x = Xi — (-D r V r )| x=Ii ) 

with a similar expression in the right domain. The differentiation matrices are accurate to the design order 
of the method. Thus, the first three terms to the right of the equality, reduce to the truncation error of the 
spatial approximation. (Except for a finite number of points that are lower by 1 order near the interfaces 
and the physical boundary). Examining the truncation error from the penalty terms, we observe that V 
is smooth across the interface, and Vj, - V ri = 0. Thus, we only need that \ x=Xi and D r V rx=x . 
approximate the first derivative to the design order of accuracy. The exact nature of the solution error near 
the boundaries is extremely complicated due to the points treated less accurately in that vicinity. More 
details will be presented in a future work on this subject. We show by numerical example, however, that 
order reduction occurs when the interface derivative is treated with less than design accuracy. (See Table 5). 


5.1. Uniform Grid.. Now we demonstrate that the physical boundary conditions must be imposed 
with accuracy of at least design order to maintain global design accuracy. This condition is a natural 
consequence of the overall dependence of the solution on the boundary conditions. The test problem we use 
is the Burgers’ equation 


(5.1) 


dU dU _ d 2 U 
dt + dx dx 2 


- 1 < x < 1 , t > 0, 


with the exact solution 


(5.2) U(x,t) = -atanh(a X °* ) 4 c, -oo < x < oo,t < 0. 

The solution of (5.1) requires imposition of boundary conditions at each end of the physical domain. We 
choose Robin boundary conditions of the form 


du , 


au(-l y t) - = 9 




du 

7 u(l,t) - <5 — |i = gi(t) 


such that the problem is mathematically well-posed. (See Hesthaven and Gottlieb [2] for the constraints 
on a, /?, 7 , and £). The physical boundary conditions were imposed in penalty form, as described in the 
work of Hesthaven and Gottlieb [2]. The time- advancement scheme is the five-stage fourth-order low-storage 
Runge-Kutta scheme. The time step was chosen to ensure that the temporal error in the formulation was 


8 



small relative to the spatial error. The simulation is run to a physical time of T — 1, and the viscosity is 
determined by the value e = 5 1CT 1 . 

Tables 1 to 4 show the results of a grid-refinement study on a single domain with a fourth-order explicit 
interior scheme. The accuracy of the boundary closure and of the physical boundary condition are parameters 
in the study. Table 2 shows the results of the refinement study with a uniformly fourth-order-accurate scheme 
(4, 4-4-4, 4) with the derivative term in the Robins’ boundary conditions approximated to 0(Ax 4 ). We note 

Table 5.1 

L 2 Solution Errors: Convergence rate of uniformly fourth- order sch eme 


N 

LOG terror 

Rate 

33 

-3.847 


65 

-4.082 

2.31 

129 

-5.239 

3.84 

257 

-6.486 

4.14 

513 

-7.731 

4.14 

1025 

-8.960 

4.87 


that the convergence rate in Table 1 is fourth order and that the design accuracy is achieved. 

Table 2 shows the second study in which boundary closure accuracy is relaxed by one order. The 
resulting scheme (3, 3-4-3, 3) is third order locally at each boundary and fourth order in the interior. (Both 
the inviscid and viscous stencils are reduced by one order of accuracy near the boundaries.) The physical 
boundary condition is still approximated to 0( Ax 4 ). 

Table 5.2 

L 2 Solution Errors: Convergence rate of fourth- order scheme with third-order closure at boundaries. 


N 

LOG terror 

Rate 

33 

-3.694 


65 

-4.797 

3.66 

129 

-5.971 

3.90 

257 

-6.117 

3.81 

513 

-7.276 

3.85 

1025 

-9.455 

3.92 


We note that the convergence rate in Table 2 asymptotes to fourth order and that the absolute levels of 
error are comparable to those obtained using the (4, 4-4-4, 4) scheme. Again, design accuracy is achieved. 

Table 3 shows the third study, in which boundary closure accuracy is relaxed by two orders. The resulting 
scheme (2, 2-4-2, 2) is second order locally at each boundary and fourth order in the interior. (Only the viscous 
terms are reduced by two orders of accuracy near the boundaries.) The physical boundary condition is still 
approximated to 0( Aa: 4 ). 

We note that the convergence rate in Table 3 asymptotes to third order, which is a reduction in global 
accuracy of one order. This behavior is consistent with Gustafsson’s [13] theory, specifically, that global 
solution accuracy allows a finite number of stencils to be reduced by one order of accuracy. 

Table 4 shows the final study, in which boundary closure accuracy is uniformly fourth-order accurate 
(4, 4-4-4, 4). The physical boundary condition is approximated to 0{ Ax 3 ), however. The convergence rate in 



Table 5.3 

Z/2 Solution Errors: Convergence rate of fourth-order scheme with second-order closure at boundaries. 


N 

LOGioerror 

Rate 

33 

-2.974 


65 

-4.074 

3.65 

129 

-5.519 

4.80 

257 

-6.284 

2.54 

513 

-7.048 

2.54 

1025 

-7.898 

2.82 


Tabic 4 asymptotes to third order, which is a reduction in global accuracy by one order. 

Table 5.4 

L 2 Solution Errors: Convergence rate of uniformly fourth-order scheme, using third-order accurate boundary conditions. 


N 

LOGioerror 

Rate 

33 

-3.004 


65 

-4.002 

3.32 

129 

-4.764 

2.53 

257 

-5.636 

2.90 

513 

-6.531 

2.97 

1025 

-7.898 

2.82 


This series of tests on the single domain indicates the need to impose the physical boundary condition 
with design accuracy. However, closing the near boundary stencils with an accuracy that is one order 
less than the design interior accuracy appears to be sufficient. A similar conclusion was reached with a 
sccond-ordcr-accurate scheme (1-2-1) and second-order physical boundary conditions. 

We now demonstrate by numerical example that these results generalize to the case of multiple domains. 
Table 5 shows a grid-refinement study that compares one and eight spatial domains. The numerical test 
problem is the previously described Burgers’ equation using a value of e = 10 2 . The numerical scheme used 
in both cases is the (3, 3, 3, 3-4-3, 3, 3, 3) scheme with physical boundary conditions imposed to an accuracy of 
0( Ax 4 ). 

Table 5.5 

L 2 Solution Errors: Convergence rate of fourth-order scheme with third-order closure at interfaces, on multiple domain 
problem. 



1 domains 


8 domains 


N 

LOGioerror 

Rate 

LOGioerror 

Rate 

97 

-2.148 


-2.125 


193 

-3.016 

2.88 

-3.143 

3.38 

385 

-4.214 

3.98 

-4.485 

4.45 

769 

-5.372 

3.85 

-5.656 

3.38 

1537 

-6.505 

3.76 

-6.866 

4.02 

3063 

-7.664 

3.85 

-8.055 

3.95 
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We note that the convergence rate in Table 5 asymptotes to fourth order, for both the one- and eight- 
domain cases. This example demonstrates that design accuracy is achieved with multiple domains so long 
as the physical boundary conditions are imposed with design accuracy and the numerical closures near the 
interfaces are at most one order of accuracy less than the design accuracy of the interior scheme. 
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Fig. 5 . 1 . The Burgers equation solved using a sixth-order scheme with randomly generated interface points. 


5.2. Nonuniform Domain. The final problem we solve is the nonlinear Burgers’ equation with un- 
equally spaced subdomains and a sixth-order scheme. Details of the numerical discretization arc included in 
the appendix. The Burgers’ equation in the form of equation (5.1) is solved throughout the domain with a 
viscosity parameter of c = 10 -2 . The domain is divided into 12 subdomains, each with the same number 
of points and a uniform local discretization. The domain interfaces are placed randomly throughout the 
domain. The ratio of maximum to minimum subdomain size is about 15:1. Figure 1 shows the solution at 
four different times. The “symbols” at the top of the figure show the positions of the 11 interface points. 
The profiles are smooth and monotone for this discretization. Figure 2 shows the logarithm of the solution 
error plotted as a function of space on the sequence of five grids. 

This problem demonstrates the stability and accuracy of the new interface treatments. The discretiza- 
tions asymptote to a convergence rate of sixth order on the sequence of grids. Table 6 shows the convergence 
rate of the calculations, for two different values of the parameter e . The steep gradients are resolved to 
high-order on all grids for e = 10 -2 . For e = 2 10~ 3 , the two coarsest grids are not yet achieving high-order 
accuracy, and two-point grid oscillations exist in the solution. Further reduction of e causes numerical insta- 
bility, emanating from the interface location, as the gradients pass the interface. Increasing the robustness 
of the interface conditions for marginally resolved/discontinuous cases is the focus of current research. 

6. Conclusions. We focus on high-order finite difference schemes, which satisfy the summation-by- 
parts (SBP) discretization framework. We show stable and conservative interface treatments of arbitrary 
spatial accuracy for the linear ad vect ion- diffusion equation. Problems with multiple domains and abruptly 
changing mesh sizes are considered. 
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Fig. 5.2. Errors obtained from Burgers equation solved on a sequence of grids with a sixth- order scheme. 

Table 5.6 

L 2 Solution Errors: Convergence of sixth-order scheme with twelve subdomains and interfaces distributed randomly. 


N 

e = nr 2 

LOGioerror 

Rate 

£ = 21(T 3 
LOGioerror 

Rate 

145 

-3.090 


-1.376 


289 

-4.641 

5.15 

-1.865 

1.62 

577 

-5.915 

4.22 

-3.053 

3.95 

1153 

-7.520 

5.33 

-4.574 

5.05 

2305 

-9.370 

6.15 

-5.834 

4.18 


Finite- difference operators are shown to admit design accuracy {p th - order global accuracy) when p — I th - 
order stencil closures are used near boundaries if the physical boundary conditions are imposed with p tAl -order 
accuracy. Finite-difference operators of up to sixth order arc constructed which satisfy the constraints of the 
new interface procedures. 

Accurate sixth order calculations are achieved for the nonlinear Burgers equation on a twelve subdomain 
problem having randomly distributed interfaces. 

REFERENCES 

[1] D. A. Kopriva, Spectral methods for the Euler equations-The blunt body problem revisited , AIAA 

Journal, 29, (1991), pp. 1458-1462. 

[2] J. S. Hesthaven and D. Gottlieb, A Stable Penalty Method for the Compressible Naxner- Stokes 

Equations : I. Open Boundary Conditions , SIAM J. Sci. Comput. 17, 3, (1996), pp 579-612. 

[3] H.-O. Kreiss and G. Scherer, Finite element and finite difference methods for hyperbolic partial 

differential equations , Mathematical Aspects of Finite Elements in Partial Differential Equations, 


12 



Academic Press, New York, (1974). 

[4] P. OLSSON, High-order difference methods and data-parallel implementation , PhD Thesis, Uppsala 

University, Department of Scientific Computing, (1992). 

[5] P. OLSSON, Summation by Parts , Projections, and Stability /, Math. Comp., 64, (1995) pp. 1035-1065. 

[6] P. OLSSON, Summation by Parts , Projections , and Stability //, Math. Comp., 64, (1995) pp. 1473-1493. 

[7] B. STRAND, High- Order Difference Approximations for Hyperbolic Initial Boundary Value Problems , 

PhD Thesis, Uppsala University, Department of Scientific Computing, (1996). 

[8] B. Strand, Summation by Parts for Finite Difference Approximations for d/dx , J. Comp. Phy., 110, 

No. 1, (1994), pp. 47-67. 

[9] M. H. Carpenter, D. Gottlieb, and S. Abarbanel, The Stability of Numerical Boundary Treat- 

ments for Compact High-Order Finite- Difference Schemes , J. Comp. Phy., 108, No. 2, (1993), pp. 
272-295. 

[10] M. H. Carpenter and D. Gottlieb, Spectral Methods on Arbitrary Grids , J. Comp. Phy., 129, No. 

0234, (1996), pp. 74-86. 

[11] M. H. Carpenter and J. Otto, High-Order Cyclo-Difference Techniques: An Alternative to Finite 

Differences J. Comp. Phy., 118, (1995), pp. 242-260. 

[12] P. D. Lax and B. Wendroff, Systems of Conservation Laws , Comm. Pure Appl. Math, 13, (1960), 

pp. 217-237. 

[13] B. GUSTAFSSON, The Convergence Rate for Difference Approximations to Mixed Initial Boundary Value 

Problems , Math. Comp. 29, 130, (1975), pp. 396-406. 

Appendix I. Stability. Here we show the algebra involved in proving THEOREM 3.1. We begin by 
restating of the stability condition presented in equation (3.6), governing the total energy of the system: 

(- 1 ) + IMIpJ ^ 

where w; = [ Ui,Vi , (Diu) i , (_D r n)J, and the boundary matrix defined in equation (3.7) is defined by 

-(<ti+<73) c(1+<7 2 ) -ea 2 
a 4- 2 <T 3 — 604 e( — 1 + 04 ) 

—6(74 —2 eai 0 

e( — 1+0-4) 0 — 2ea r 

The stability of this matrix is easier to analyze if it is rotated with a similarity transformation. Define 
the new vector w = S w such that: 

U{ — 1—10 0 

Ui + Vi 1110 0 

(D l u) i -(D r v) l “72 0 0 1 -1 

{D t u) i + (D r v) J [ooil 

The similarity rotation matrix has the property S T = .S' ~ 1 as can easily be verified. The rotation matrix 
5 can be used to transform the stability condition defined by equation (3.6) into the following equivalent 
condition: 

(.4) wfAPwi = w'[S T SM l S T Sw l = w T M l w<0; 


Ui 

Vi 

( Diu\ 
{D r v\ 
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where 


(.5) 


2(<7i 4 < 73) 
-(-0-1 + <r 3 4 a) 
e(a 2 + cr 4 ) 


e 


( <7\ + cr 3 + a) 
0 

-c(-a 2 4- <J 4 ~ 1) 
0 


e(<7 2 -f <r 4 ) 
-c(-<r 2 + cr 4 - 1) 
-e(a r +a/) 
e(a r - cq) 


c 

0 


e(a r - a*) 
-e(a r 4 a/) 


To ensure negative definiteness, every sub-matrix in the matrix M l must be negative definite. We 
observe by inspection that (cr\ + (7 3 ) < 0 is a necessary condition. Analyzing the 2x2 sub-matrices along the 
diagonal, we obtain the necessary conditions ( — cr 4 -f a 3 4 a) = 0, and c(— cr 2 4 cr 4 — 1) = 0. Substituting 
the equalities (— o\ 4 <j 3 4 a) = 0 and (— <7 2 4 rr 4 — 1) = 0 into the matrix M x yields: 


(. 6 ) 


M l = 


2(2 cti — a) 
0 

c(2<j 2 4 1) 


0 

0 

0 

0 


e(2cr 2 + 1) 
0 

-e(a r 4 af) 
e(a r - oq) 


e 

0 


e(a r - a t ) 
-e(a r +at) 


A symmetric matrix can be rotated into diagonal form by an orthogonal matrix, making the condition 
of negative semi-definiteness 


w T U T D l Uvr < 0; 

where U is the orthogonal matrix that satisfies U T D l U ~ M\ Pre- and post- multiplication of M* by 
suitable rotation matrices M\ = R\ 7 M X R\, yield the equivalent condition 

w T R 1 T U T D i UR l w < 0; 

The matrix R \ , chosen to yield a diagonal expression for the matrix M\ is 

0 0 0 

1 0 0 

0 1 0 

0 £4,3 1 

with 

r _ -c(2tT2 4 1) 

34 “ 2(2(7 i-a) 

L _ — 2e(a r a 2 - aia 2 4a r ) 

44 e(4crf 4 4<7 2 4 1)4 (4<Ji - 2a)(a r 4 a/) 

L — 4“ 1) 4 (— 4cri 4 2a) (a r — ai)) 

4,3 €(4a| 4 4 <t 2 4 1)4 (4<7i - 2a) (a r 4 a/) 

The diagonal elements of M\ are 


(.7) 


L v 


1 

0 

£3,1 

T44 


Ai = 2(2(7 1 — a) 
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A2 = 0 

= ~ f ( £ ( 4 < 7 1 + 402 + 1 ) + ( 4 < 7 x - 2 a) (ay + g t )) 
2(2<7i — a) 

_ -4e(a r e(g 2 + l) 2 + cqecrf + (4 eg - 2a)cqg r ) 
e (4 (T 2 4" 4(72 + 1) + (4(71 ~ 2a) (a r + a/) 


These eigenvalues must be less than or equal to zero to ensure stability of the interface condition. The 
resulting condition of stability is 


Combining this expression with the constraints <r 3 = <n -a and a 4 = <r 2 +l yield the conditions of THEOREM 

Appendix II. Stencils. Wc now present the specific form of the stencils that satisfy the SBP stability 
requirements, and the accuracy requirements shown necessary in the previous numerical study. At second 
order, the discretization matrix for the advcction terms that satisfy the constraint Aj = P~ l Q is 


(. 8 ) 


where 


-2 2 

-1 0 


A = 


1 

2Ax 


-1 0 1 
-2 2 




" 1 
2 

1 


-1 1 

-1 0 1 


(.9) 

P= Ax 

1 

1 

2 _ 

; Q=l 


-1 0 1 
— 1 i 


The discretization matrix for the diffusion terms that satisfies the constraint A 2 = P~ l (-S T R + D)S 

is ' 


(. 10 ) 


1 -2 1 

1 -2 1 



1 -2 1 

1 -2 1 
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where 


(.id 


and 


(. 12 ) 



The matrix R can be shown to be positive definite (and symmetric). 

The fourth-order discretization that satisfies the SBP constraints was originally derived in the work of 
Strand [8]. The coefficients rl and r2 below are different from those proposed by Strand and are chosen so 
that the resulting discretization + = P l Q has the standard four-point third-order stencil at the first grid 
point. The values of rl and r2 are 


(.13) 


-(21 77^295369 - 1166427) 
25488 

(66195 v ^v / j^ - 35909375) 
101952 


and the matrices P and Q are 


P — Ax 


' —(216 r2+2160 rl— 2125) 
12960 

(81 r2+675 rl+415) 
540 

-(72 r2+720 r 1+445) 
1440 

-(108 r2+756rl+421) 
1296 


(81 r2+675 rl+415) 

540 

-(4104 r 2+ 32400 r 1 + 1 1225) 


4320 

(1836 r2+ 14580 rl + 7295) 
2160 

-(216 r2+2160 r 1+655) 


4320 


-(72 r2+72Q rl+445) 
1440 

(1836 r2+ 14580 rl + 7295) 
2160 

-(4104 r 2+ 32400 rl+12785) 
4320 

(81 r2+675 rl+335) 

540 


-(108 r2+756 rl + 421) 

1296 

-(216 r2+2160 rl+655) 

4320 

(81 r2+675 rl+335) 

540 

-(216 r2+2160 rl-12085) 
12960 

1 


(. 14 ) 
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and 


Q= 




2 

(864 r2+ 6480 rl+305) 
4320 

— (216 r2-f- 1620 r l-j-725) 


540 

(864 r2+6480 rl-f-3335) 
4320 


-(864 r2+6480 rl+305) 
4320 

0 

(864 r2+6480 rl + 2315) 
1440 

-(108r2+810rl+415) 

270 


(216r2+1620rl + 725) 

-(864 r2-f 6480 rl+3335) 



540 

-(864 r2+6480 rl + 2315) 

4320 

(108r2+810rl+415) 



1440 

0 

(864 r2+6480 rl+785) 
4320 
1 

12 

270 

-(864 r2+6480 rl+785) 
4320 

0 

-8 

12 

-1 

12 

8 

12 

0 

-1 

12 

8 -1 

12 12 


(.15) 


Only the inflow boundary portion of the matrices P and Q is shown. The outflow coefficients are the negative 
transpose of the inflow coefficients. The matrix P is symmetric and positive definite. 

The discretization matrix for the diffusion terms that satisfies the constraint A 2 = P~ 1 (-S T R + D)S 

is: ' 


(.16) 


where 


(.17) 


35 —26 19 —14 ii 

12 3 2 ~ 12 

11 =5 1 1-1 

12 3 2 3 72 

1 16 -30 16 -1 

12 12 12 12 72 


A = 


1 

(Axf 


-1 

16 

-30 

16 

-1 

12 

12 

12 

12 

12 

-1 

1 

1 

-5 

11 

12 

3 

2 

3 

12 

11 

-14 

19 

-26 

35 

12 

3 

2 

3 
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The matrix R is too complicated to report here but can be shown to be positive definite. This numerical 

scheme is referred to as (3,3, 3, 3-4-3, 3,3, 3), which denotes the fact that the four points nearest to the boundary 
are closed with third-order formulas. 

The sixth-order discretization that satisfies the SBP constraints was originally derived in the work of 
Strand [8], The coefficients rl, r2, and r3 below arc different from those proposed by Strand and are chosen 
so that the resulting discretization A 1 = P^Q has the standard six-point fifth-order stencil at the first grid 
point. This choice produces remarkably good stability characteristics at the boundary. The coefficients are 


(. 18 ) 


rl = -3.6224891259957 
r2 = 96.301901955532 
r3 = -609.5813881563 
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The symmetric P and nearly skew- symmetric Q matrices have the entries A x - P l Q, where 


(.19) 


p(M) = 
P( 1.2) = 
P( 1.3) = 

P(M) = 

p(l, 5 ) = 

P(l>6) = 

P( 2,2) = 

P( 2,3) = 

P(2,4) = 

p(2,5) = 

P(2,6) = 

P(3,3) 

P(3,4) 

p(3,5) 

p(3,6) 

p(4,4) 

p(4,5) 

p(4,6) 

P(5,5) 

P( 5,6) 

p(6,6) 

g(i, i) 
9 ( 1 , 2 ) 
9(1,3) 
9(1,4) 


-(14400 r2 + 30 2400 rl -7420003) 

36288000 

-(75600 r3 + 1497600 r2 + 1 1944800 rl - 59330023) 
21722800 

-(9450 r3 + 202050 r2 + 1776600 rl - 7225847) 

340200 

(900 r2 + 18900 rl -649) 

226800 

(86400 r3 + 1828800 r2 + 15854400 rl - 6615 0023) 
3110400 

(378000 r3 + 7747200 r2 + 65167 200 rl - 279318239) 
188640000 

(302400 r3 + 6091200 r2 + 49896000 rl - 210294289) 
7257600 

(3780 r3 + 82575 r2 + 741 825 rl - 2991977) 

34020 

(5400 r3 + 104400 r2 + 810000 rl - 3756643) 

" 129600 

-(529200 r3+ 11107200 r2 + 95508000 rl -400851749) 
2419200 

(86400 r3 + 1828800 r2 + 15854400 rl - 65966279) 
3110400 

- (51300 r3 + 1094400 r2 + 9585000 rl - 39593423) 
64800 

(120960 r3 + 2584800 r2 + 22680000 r l - 93310367) 
181440 

(5400 r3 + 104400 r2 + 810000 rl - 3766003) 

129600 

(900 r2 + 18900 rl - 37217) 

226800 

-(17100 r3 + 364800 r2 + 3195000 rl - 13184701) 
21600 

(3780 r3 + 82575 r2 + 741825 rl - 2976857) 

34020 

-(1890 r3 + 40410 r2 + 355320 rl - 1458223) 

: 68040 

(302400 r3 + 6091200 r2 + 49896000 rl - 213056209) 

: 7257600 

-(75600 r3 + 1497600 r2 + 1194480 0rl - 54185191) 

: 21722800 

-(14400 r2 + 302400 rl - 36797603) 

= 36288000 

- (- 1 ) 

2 

(415800 r3 + 8604000 r2 + 72954000 rl - 283104553) 

= 32659200 

(120960 r3 + 2672640 r2 + 241 92000 rl - 100358119) 

= 6531840 

-(25200 r3 + 542400 r2 + 478 8000 rl - 19717139) 

= ~ 403200 
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9 ( 1 , 5 ) 

9(1,6) 

9 ( 2 , 2 ) 

9 ( 2 , 3 ) 

9 ( 2 , 4 ) 

9 ( 2 , 5 ) 

9(2,6) 

9 ( 3 , 3 ) 

9 ( 3 , 4 ) 

9 ( 3 , 5 ) 

9 ( 3 , 6 ) 

9 ( 4 , 4 ) 

9 ( 4 , 5 ) 

9 ( 4 , 6 ) 

9 ( 5 , 5 ) 

9 ( 5 , 6 ) 

9 ( 6 , 6 ) 


(604800 r3 + 13363200 r2 + 120960000 rl - 485628701) 
32659200 

(41580 r3 + 860400 r2 + 7295400 rl - 31023481) 

3265920 

0 

- (9450000 r3 + 200635200 r2 + 17471 16000 rl - 7286801279) 

_____ 

(21 168000 r3 + 449049600 r2 + 3907008000 rl - 16231108387) 

32659200 

-(165375 r3 + 3516300 r2 + 30665250 rl - 126996371) 

453600 

(604800 r3 + 13363200 r2 + 120960000 rl - 482536157) 
32659200 

0 

-(6993000 r3 + 148096800 r2 + 1286334000 rl - 5353075351) 

______ 

(21168000 r3 + 449049600r2 + 3907008000rl - 16212561187) 

32659200 

-(75600 r3 + 1627200 r2 + 14364000 rl - 58713721) 

1209600 

0 

- (9450000 r3 + 200635200 r2 + 17471 16000 rl - 7263657599) 

32659200 

(604800 r3 + 13363200 r2 + 120960000 rl - 485920643) 
32659200 

0 

(415800 r3 + 8604000 r2 + 72954000 rl - 286439017) 


The matrix P is symmetric and positive definite for this choice of parameters. 

The discretization matrix for the diffusion terms that satisfies the constraint An — P~ x 
is 


(. 20 ) 


where 


(. 21 ) 


A = 


1 

180 (Ax) 2 


Q— 1 
Ax 


1 iM 
20 


+812 

—3132 

+5265 

-5080 

+2970 

-972 

+ 137 

+ 137 

-147 

-255 

+470 

-285 

+93 

-13 

-13 

+228 

-420 

+200 

+15 

-12 

+2 

2 

-27 

270 

-490 

270 

-27 

2 


(- 15 ) 20 

(- 15 ) 6 (- 1 ) 1 


"-1 

2 3 

4 5 6 




, D = 

0 






The matrix R is too complicated to report here but can be shown to be 


positive definite. 


(~S r R + D)S 
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